Fully reported case
Simulation
params_b_rw_FR <- list(
data = list(
alpha_lamb = alpha_lamb,
beta_lamb = beta_lamb,
T = T,
date_start = as.Date("2024-01-01"),
D = D
),
q_model = list(
method = "b_rw",
method_params = list(mu_logb = log(0.7), sigma_logb = 0.1, mu_logitphi = 1, sigma_logitphi = 0.1)
)
)
b_rw_FR <- simulateData(params_b_rw_FR)
par(mfrow = c(2, 1))
plot(b_rw_FR$b, pch = 19, type = "b")
plot(b_rw_FR$phi, pch = 19, type = "b")
par(mfrow = c(1, 1))
matplot(t(b_rw_FR$q), type = "l", lty = 1, ylim = c(0, 1))
Exploratory analysis
# exploritary analysis
page_num <- ceiling(nrow(b_rw_FR$case_reported_cumulated)/16)
exp_plot_b_rw <- fit_exp_plot(b_rw_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# print(exp_plot_b_rw)
#
# exp_plot_b_rw$coefficients
exp_b_data_b_rw<- data.frame( date = as.Date(rownames(b_rw_FR$case_reported_cumulated)),
b = exp_plot_b_rw$coefficients$b)
exp_b_plot_b_rw <- ggplot(exp_b_data_b_rw, aes(x = date, y = b)) +
geom_point(color = "black", size = 1.5) +
geom_smooth(method = "loess", se = TRUE,
color = "blue", fill = "grey", alpha = 0.5) +
theme_minimal() +
labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")
print(exp_b_plot_b_rw)
## `geom_smooth()` using formula = 'y ~ x'
Model fitting
# out_b_rw_FR <-
# nowcasting_moving_window(b_rw_FR$case_reported_cumulated, scoreRange = scoreRange,
# case_true = b_rw_FR$case_true,
# start_date = as.Date("2024-01-01"),
# D = D,
# methods =models_to_use,
# compiled_models = compiled_models,
# iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
# num_chains = 3, suppress_output = T,
# posterior_draws_path = file.path(posterior_draws_path, "b_rw")
# )
#
# save(out_b_rw_FR, file = file.path(data_save_path, "FR_b_rw.RData"))
load( file.path(data_save_path, "FR_b_rw.RData"))
results_b_rw_FR <- nowcasts_table(out_b_rw_FR, D = D, report_unit = "day",
methods = models_to_use)
plots_b_rw_FR <- nowcasts_plot(results_b_rw_FR, D = D, report_unit = "day", methods = models_to_use)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
for (i in 1:length(results_b_rw_FR)) {
print(calculate_metrics(results_b_rw_FR[[i]],
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 77.76 62.00 53.53 60.54 120.61 0.6 q_constant
## 2 15.48 8.05 8.07 6.66 57.31 1.0 b_constant
## 3 16.19 8.09 8.20 6.71 67.21 1.0 b_rw
## 4 16.11 8.18 8.19 6.75 61.61 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 67.30 12.43 46.87 11.57 107.16 0.75 q_constant
## 2 31.90 4.50 15.75 3.07 80.65 1.00 b_constant
## 3 27.66 3.88 12.12 2.40 95.06 1.00 b_rw
## 4 32.36 4.49 15.30 2.93 90.51 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 54.52 4.65 34.19 3.75 109.34 0.87 q_constant
## 2 30.63 2.88 15.09 2.07 103.67 0.97 b_constant
## 3 22.12 1.91 8.14 1.12 124.34 1.00 b_rw
## 4 22.88 2.04 9.28 1.27 124.95 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 24.40 2.43 15.79 1.82 114.33 1 q_constant
## 2 17.62 2.33 9.94 1.59 112.18 1 b_constant
## 3 19.45 2.14 7.97 1.10 126.29 1 b_rw
## 4 13.36 1.71 6.42 1.02 125.13 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 33.60 6.90 15.32 2.61 110.07 0.98 q_constant
## 2 30.66 6.63 11.12 2.32 108.13 0.98 b_constant
## 3 15.13 3.31 5.80 1.24 112.65 1.00 b_rw
## 4 19.68 4.35 6.14 1.40 113.17 1.00 b_ou
for (i in 1:length(results_b_rw_FR)) {
print(calculate_metrics(data.table::last(results_b_rw_FR[[i]],D_check),
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 109.38 67.98 96.61 66.08 211.23 0.2 q_constant
## 2 21.86 9.17 15.36 8.42 97.01 1.0 b_constant
## 3 22.86 9.31 15.58 8.40 116.41 1.0 b_rw
## 4 22.75 9.33 15.60 8.43 105.81 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 117.55 15.88 109.96 14.80 226.63 0.2 q_constant
## 2 61.62 7.75 49.53 6.40 166.02 1.0 b_constant
## 3 53.39 6.65 37.18 4.76 219.22 1.0 b_rw
## 4 62.92 7.89 49.07 6.33 202.02 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 110.77 8.05 102.23 7.52 231.61 0.4 q_constant
## 2 68.07 4.85 58.24 4.25 220.60 0.8 b_constant
## 3 52.28 3.67 34.30 2.43 337.63 1.0 b_rw
## 4 54.72 3.83 42.74 3.04 339.24 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 46.65 4.58 35.68 3.39 198.42 1 q_constant
## 2 38.92 3.77 28.06 2.64 192.42 1 b_constant
## 3 53.08 5.24 39.78 3.83 293.43 1 b_rw
## 4 35.47 3.46 28.97 2.77 289.23 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 97.61 21.08 62.40 13.01 140.01 0.8 q_constant
## 2 93.36 20.19 56.80 11.93 137.80 0.8 b_constant
## 3 45.81 9.84 31.14 6.40 169.40 1.0 b_rw
## 4 61.00 13.19 38.42 7.95 175.02 1.0 b_ou
print(plots_b_rw_FR)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Parameter Checking
# try the third one, "2024-01-30"
T_test = T * 3/6
# q_constant
varnames_b_rw <- out_b_rw_FR$b_rw[[3]]$summary()$variable
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("sigma_log_b"), prob_outer = 0.95)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("sigma_logit_phi"), prob_outer = 0.95)
param_true = tibble(
parameter = grep("^b\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_FR$b[1:T_test]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("b"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^phi\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_FR$phi[1:T_test]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("phi"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_FR$lambda[1:T_test]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_FR$q[10,]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE)), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_rw_FR$b_rw[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
Simulation
params_b_rw_NFR <- list(
data = list(
alpha_lamb = alpha_lamb,
beta_lamb = beta_lamb,
T = T,
date_start = as.Date("2024-01-01"),
D = D
),
q_model = list(
method = "b_rw",
method_params = list(mu_logb = log(0.1), sigma_logb = 0.1, mu_logitphi = 1.5,
sigma_logitphi = 0.1)
)
)
b_rw_NFR <- simulateData(params_b_rw_NFR)
par(mfrow = c(2, 1))
plot(b_rw_NFR$b, pch = 19, type = "b")
plot(b_rw_NFR$phi, pch = 19, type = "b")
par(mfrow = c(1, 1))
matplot(t(b_rw_NFR$q), type = "l", lty = 1, ylim = c(0, 1))
b_rw_NFR$case_true - b_rw_NFR$case_reported_cumulated[,16]
## [,1]
## 2024-01-01 1
## 2024-01-02 1
## 2024-01-03 3
## 2024-01-04 1
## 2024-01-05 3
## 2024-01-06 1
## 2024-01-07 13
## 2024-01-08 11
## 2024-01-09 17
## 2024-01-10 28
## 2024-01-11 36
## 2024-01-12 41
## 2024-01-13 41
## 2024-01-14 43
## 2024-01-15 41
## 2024-01-16 50
## 2024-01-17 41
## 2024-01-18 57
## 2024-01-19 57
## 2024-01-20 59
## 2024-01-21 69
## 2024-01-22 79
## 2024-01-23 100
## 2024-01-24 145
## 2024-01-25 131
## 2024-01-26 194
## 2024-01-27 217
## 2024-01-28 269
## 2024-01-29 250
## 2024-01-30 165
## 2024-01-31 235
## 2024-02-01 192
## 2024-02-02 222
## 2024-02-03 357
## 2024-02-04 377
## 2024-02-05 339
## 2024-02-06 380
## 2024-02-07 304
## 2024-02-08 348
## 2024-02-09 378
## 2024-02-10 281
## 2024-02-11 281
## 2024-02-12 260
## 2024-02-13 268
## 2024-02-14 263
## 2024-02-15 184
## 2024-02-16 168
## 2024-02-17 143
## 2024-02-18 133
## 2024-02-19 116
## 2024-02-20 106
## 2024-02-21 95
## 2024-02-22 65
## 2024-02-23 61
## 2024-02-24 25
## 2024-02-25 6
## 2024-02-26 1
## 2024-02-27 4
## 2024-02-28 10
## 2024-02-29 4
Exploratory analysis
# exploritary analysis
page_num <- ceiling(nrow(b_rw_NFR$case_reported_cumulated)/16)
exp_plot_b_rw <- fit_exp_plot(b_rw_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
# print(exp_plot_b_rw)
#
# exp_plot_b_rw$coefficients
exp_b_data_b_rw<- data.frame( date = as.Date(rownames(b_rw_NFR$case_reported_cumulated)),
b = exp_plot_b_rw$coefficients$b)
exp_b_plot_b_rw <- ggplot(exp_b_data_b_rw, aes(x = date, y = b)) +
geom_point(color = "black", size = 1.5) +
geom_smooth(method = "loess", se = TRUE,
color = "blue", fill = "grey", alpha = 0.5) +
theme_minimal() +
labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")
print(exp_b_plot_b_rw)
## `geom_smooth()` using formula = 'y ~ x'
Model fitting
# out_b_rw_NFR <-
# nowcasting_moving_window(b_rw_NFR$case_reported_cumulated, scoreRange = scoreRange,
# case_true = b_rw_NFR$case_true,
# start_date = as.Date("2024-01-01"),
# D = D,
# methods =models_to_use,
# compiled_models = compiled_models,
# iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
# num_chains = 3, suppress_output = T,
# posterior_draws_path = file.path(posterior_draws_path, "b_rw")
# )
#
# save(out_b_rw_NFR, file = file.path(data_save_path, "NFR_b_rw.RData"))
load( file.path(data_save_path, "NFR_b_rw.RData"))
results_b_rw_NFR <- nowcasts_table(out_b_rw_NFR, D = D, report_unit = "day",
methods = models_to_use)
plots_b_rw_NFR <- nowcasts_plot(results_b_rw_NFR, D = D, report_unit = "day", methods = models_to_use)
for (i in 1:length(results_b_rw_NFR)) {
print(calculate_metrics(results_b_rw_NFR[[i]],
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 25.18 13.66 13.34 11.68 83.11 1.0 q_constant
## 2 81.13 50.91 53.68 50.29 34.20 0.1 b_constant
## 3 81.50 50.76 53.87 50.10 35.71 0.1 b_rw
## 4 80.31 50.53 53.15 49.92 35.60 0.1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 119.63 19.24 75.22 16.97 92.81 0.60 q_constant
## 2 129.31 21.16 83.42 19.22 82.61 0.55 b_constant
## 3 66.77 11.81 38.68 9.65 175.16 1.00 b_rw
## 4 76.49 13.18 40.90 10.63 170.07 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 183.66 17.40 121.32 15.64 101.94 0.40 q_constant
## 2 122.54 10.38 65.11 7.48 112.74 0.83 b_constant
## 3 176.11 14.76 85.15 10.05 177.45 0.87 b_rw
## 4 187.74 15.32 89.54 10.13 160.38 0.80 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 249.41 21.24 180.41 18.36 103.21 0.30 q_constant
## 2 154.88 13.37 106.39 10.93 116.46 0.55 b_constant
## 3 85.01 9.56 49.94 6.37 231.75 0.98 b_rw
## 4 126.04 12.41 73.58 7.98 209.96 0.88 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 212.52 20.39 159.93 17.66 99.91 0.30 q_constant
## 2 108.60 12.29 87.94 11.30 112.79 0.36 b_constant
## 3 41.08 7.25 28.40 4.96 195.04 1.00 b_rw
## 4 54.01 8.99 35.38 5.77 188.19 0.98 b_ou
for (i in 1:length(results_b_rw_NFR)) {
print(calculate_metrics(data.table::last(results_b_rw_NFR[[i]],D_check),
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 35.55 15.73 24.95 13.57 143.22 1 q_constant
## 2 114.37 54.67 98.70 54.24 55.01 0 b_constant
## 3 114.90 54.79 99.12 54.35 58.01 0 b_rw
## 4 113.21 54.24 97.68 53.79 57.80 0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 227.40 28.63 219.14 28.19 183.62 0 q_constant
## 2 244.33 30.84 236.68 30.49 164.01 0 b_constant
## 3 129.16 16.04 120.63 15.40 443.04 1 b_rw
## 4 149.88 18.07 130.84 16.39 423.83 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 395.59 28.62 392.93 28.30 209.42 0.0 q_constant
## 2 284.81 20.75 276.31 20.01 238.01 0.2 b_constant
## 3 423.53 30.19 418.58 29.90 498.06 0.2 b_rw
## 4 452.63 32.16 446.30 31.82 423.05 0.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 335.88 31.72 332.02 30.88 169.41 0.0 q_constant
## 2 216.86 21.23 208.19 19.64 194.81 0.0 b_constant
## 3 189.36 19.85 168.41 16.52 575.30 0.8 b_rw
## 4 287.13 28.59 277.55 26.49 423.63 0.4 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 152.48 28.22 147.75 27.96 119.22 0.0 q_constant
## 2 81.47 14.80 75.28 14.10 138.21 0.6 b_constant
## 3 72.75 14.07 69.27 13.35 326.43 1.0 b_rw
## 4 116.70 22.80 114.97 22.30 252.42 0.8 b_ou
print(plots_b_rw_NFR)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Parameter Checking
# try the third one, "2024-01-30"
T_test = T * 3/6
# q_constant
varnames_b_rw <- out_b_rw_NFR$b_rw[[3]]$summary()$variable
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("sigma_log_b"), prob_outer = 0.95)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("sigma_logit_phi"), prob_outer = 0.95)
param_true = tibble(
parameter = grep("^b\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_NFR$b[1:T_test]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("b"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^phi\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_NFR$phi[1:T_test]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("phi"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_NFR$lambda[1:T_test]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_NFR$q[10,]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws(grep("^q\\[10,.+\\]$", varnames_b_rw, value = TRUE)), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_b_rw, value = TRUE),
x = b_rw_NFR$case_true[1:T_test, 1]
)
mcmc_areas(out_b_rw_NFR$b_rw[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
